// Copyright 2012 Sergio, Raul, Dorival. All rights reserved.
// Use of this source code is governed by a BSD-style
// license that can be found in the LICENSE file.

// STL
#include <iostream>
#include <cmath>

// GemLab
#include "gemlab.h"

using std::cout;
using std::endl;
using GemLab::Vec_t;
using GemLab::Mat_t;
using GemLab::Array;
using GemLab::String;
using GemLab::Geo2name;
using GemLab::Geo2nvc;
using GemLab::Geo2gnd;
using GemLab::Geo2natc;
using GemLab::Geo2func;

int main(int argc, char **argv) try
{
    bool with_err = false;

    for (int geo=0; geo<GemLab::Ngeo; ++geo) {
        cout << "------------------------- " << Geo2name[geo] << " -------------------------" << endl;

        // variables
        int nvc = Geo2nvc[geo];
        int gnd = Geo2gnd[geo];
        //Mat_t natc = Geo2natc[geo];
        Vec_t S   (nvc);
        Mat_t dSdR(nvc, gnd);

        // check S
        double tol  = 1e-17;
        double errS = 0.0;
        for (int n=0; n<nvc; ++n) {
            Vec_t rst(3);
            for (int i=0; i<gnd; ++i) {
                //rst[i] = natc(i,n);
                rst[i] = Geo2natc[geo][i][n];
            }
            Geo2func[geo](S, dSdR, rst(0), rst(1), rst(2), false);
            for (int m=0; m<nvc; ++m) {
                if (n == m) {
                    errS += fabs(S[m] - 1.0);
                } else {
                    errS += fabs(S[m]);
                }
            }
        }
        if (errS > tol) {
            cout << Geo2name[geo] << " failed with err = " << errS << endl;
            with_err = true;
        } else {
            cout << "OK" << endl;
        }
    }

    if (with_err) return 1;
    return 0;

    /*
        // check dSdR
        tol     = 1e-14
        h      := 1.0e-1
        S_temp := make([]float64, nvc)
        if name == "lin5" || name == "tri15" {
            tol = 1.0e-10
        }
        for n := 0; n < nvc; n++ {
            rst := []float64{0, 0, 0}
            for i := 0; i < gnd; i++ {
                rst[i] = natc[i][n]
            }
            // analytical
            Geo2func[geo](S, dSdR, rst[0], rst[1], rst[2], true)
            // numerical
            for i := 0; i < gnd; i++ {
                dSndRi, _ := num.DerivCentral(func(x float64, args ...interface{}) (Sn float64) {
                    rst_temp   := []float64{rst[0], rst[1], rst[2]}
                    rst_temp[i] = x
                    Geo2func[geo](S_temp, nil, rst_temp[0], rst_temp[1], rst_temp[2], false)
                    Sn = S_temp[n]
                    return
                }, rst[i], h)
                utl.Pfgrey2("  dS%ddR%d @ [% 4.1f % 4.1f % 4.1f] = %v (num: %v)\n", n, i, rst[0],rst[1],rst[2], dSdR[n][i], dSndRi)
                tol2 := tol
                if name == "tri15" && n == 11 && i == 1 {
                    tol2 = 1.0e-9
                }
                if math.Abs(dSdR[n][i] - dSndRi) > tol2 {
                    utl.Panic("%s dS%ddR%d failed with err = %g", name, n, i, math.Abs(dSdR[n][i] - dSndRi))
                }
                //utl.CheckScalar(tst, fmt.Sprintf("dS%ddR%d", n, i), tol2, dSdR[n][i], dSndRi)
            }
        }
        utl.PfGreen("OK\n")
    */
}
GEMLAB_CATCH
